
## 2.5作图合集
# 加载所需库------
if (!require("pacman")) install.packages("pacman")
pkgs <- c("dplyr", "ggplot2", "readxl", "tidyr", "lubridate", "janitor", "stringr", "writexl",
          "psych", "magrittr", "sf", "spdep", "tmap", "ggspatial","RColorBrewer","here")
pacman::p_load(pkgs, character.only = TRUE)

source("./code/2.3 人口数估发病率.R")


### 主图1：2005-2023年逐年逐月发病数------
yearly_month_incidence <- rate_month %>% 
  dplyr::select(region,year,month,incidence_rate = rate,cases = count, pop = value)

# 生成完整日期序列并填充缺失值
yearly_month_incidence1 <- yearly_month_incidence %>%
  mutate(date = ymd(paste(year, month, "01", sep = "-"))) %>%
  complete(date = seq(ymd("2005-01-01"), ymd("2023-12-01"), by = "1 month")) %>%
  replace_na(list(cases = 0))

# 步骤1：数据预处理
province_data <- yearly_month_incidence1 %>% 
  filter(region == "province") %>%  # 筛选省级数据
  mutate(date = ymd(date)) %>% 
  arrange(date)

# 计算双轴比例系数（关键）
scaling_factor <- max(province_data$cases, na.rm = TRUE) / max(province_data$incidence_rate, na.rm = TRUE)

# 可视化主题设置 -----------------------------------------------------------
idp_theme <- theme(
  text = element_text(family = "Arial", color = "#333333"),
  panel.background = element_blank(),
  panel.grid.major = element_line(color = "grey90", linewidth = 0.25),
  axis.line = element_line(color = "black", linewidth = 0.4),
  legend.position = "none",
  plot.title = element_text(size = 10, face = "bold", hjust = 0.5),
  axis.title = element_text(size = 9),
  axis.text = element_text(size = 8)
)
# 步骤2：构建双轴图
p_original <- ggplot(province_data, aes(x = date)) +
  # 主y轴：发病数柱状图
  geom_col(aes(y = cases, fill = "Cases"), 
           width = 20, 
           alpha = 0.7,
           linewidth = 0.3) +
  # 次y轴：发病率折线图
  geom_line(aes(y = incidence_rate * scaling_factor, color = "Incidence Rate"), 
            linewidth = 0.8,
            key_glyph = "path") +
  # 双轴标度设置
  scale_y_continuous(
    name = "Number of Cases",
    labels = label_number(scale_cut = cut_short_scale()),
    sec.axis = sec_axis(
      ~./scaling_factor,
      name = "Incidence\n(per 100,000)",
      labels = number_format(accuracy = 0.1)
    ),
    expand = expansion(mult = c(0, 0.05))
  ) +
  # X轴时间刻度
  scale_x_date(
    breaks = seq(ymd("2005-01-01"), ymd("2023-01-01"), by = "2 years"),
    date_labels = "%Y",
    expand = expansion(add = c(30, 60))
  ) +
  # 颜色标度（符合IDP期刊风格）
  scale_fill_manual(
    name = "",
    values = c("Cases" = "#2B8CBE"),
    labels = c("Cases"),
    guide = guide_legend(order = 1)
  ) +
  scale_color_manual(
    name = "",
    values = c("Incidence Rate" = "#D62728"),
    labels = c("Incidence Rate"),
    guide = guide_legend(order = 2)
  ) +
  # 期刊风格主题
  theme_minimal() +
  theme(
    text = element_text(family = "Arial", color = "#333333"),
    panel.grid.major = element_line(color = "grey90", linewidth = 0.25),
    panel.grid.minor = element_blank(),
    axis.line = element_line(color = "black", linewidth = 0.4),
    axis.title.y.left = element_text(size = 10,  margin = margin(r = 10)),#face = "bold",
    axis.title.y.right = element_text(size = 10,  margin = margin(l = 10)),#face = "bold",
    axis.text = element_text(size = 8),
    legend.position = "bottom",
    legend.spacing.x = unit(0.3, "cm"),
    legend.key.width = unit(1.5, "cm"),
    legend.text = element_text(size = 9),
    plot.margin = margin(15, 20, 10, 15)
  ) +
  labs(x = NULL,  title = "Original Time Series") +
  idp_theme 

# 保存高清图
# ggsave("./主图/Fig.1 2005-2023年逐年逐月发病数.png",dpi = 600)

# 加载必要包 ----------------------------------------------------------------
library(tidyverse)
library(lubridate)
library(forecast)
library(openxlsx)
library(patchwork)

# 数据预处理 ----------------------------------------------------------------
province_ts <- yearly_month_incidence1 %>% 
  filter(region == "province") %>%  # 筛选省级数据
  mutate(date = make_date(year, month, 1)) %>%  # 标准化日期
  arrange(date) %>% 
  select(date, incidence_rate = cases)

# 创建时间序列对象 ---------------------------------------------------------
ts_data <- ts(province_ts$incidence_rate,
              start = c(2005, 1),
              frequency = 12)

# STL分解 -----------------------------------------------------------------
stl_decomp <- stl(ts_data, 
                  s.window = "periodic",
                  robust = TRUE)

# 构建分解结果数据框 --------------------------------------------------------
decomp_df <- data.frame(
  date = province_ts$date,
  observed = as.numeric(stl_decomp$time.series[, "remainder"] + 
                          stl_decomp$time.series[, "trend"]),
  trend = as.numeric(stl_decomp$time.series[, "trend"]),
  seasonal = as.numeric(stl_decomp$time.series[, "seasonal"]),
  remainder = as.numeric(stl_decomp$time.series[, "remainder"])
)




# 生成各分面图 ------------------------------------------------------------
## 原始数据时序图
# p_original <- ggplot(decomp_df, aes(x = date)) +
#   geom_line(aes(y = observed), color = "#2B8CBE", linewidth = 0.6) +
#   geom_vline(xintercept = as.Date(c("2009-12-31","2015-12-01", "2019-12-01")),
#              linetype = "dashed", color = "#D62728", linewidth = 0.4) +
#   labs(x = NULL, y = "Incidence Rate (per 100,000)", 
#        title = "Original Series with Phase Division") +
#   idp_theme +
#   annotate("text", x = as.Date("2007-01-01"), y = max(decomp_df$observed)*0.8,
#            label = "Monitoring Initiation\n(2005–2009)", size = 2.8, color = "#333333") +
#   annotate("text", x = as.Date("2013-01-01"), y = max(decomp_df$observed)*0.8,
#            label = "Monitoring Enhancement\n(2010–2015)", size = 2.8, color = "#333333") +
#   annotate("text", x = as.Date("2017-03-01"), y = max(decomp_df$observed)*0.8,
#            label = "National Brucellosis Control Program\n(2016–2019)", size = 2.8, color = "#333333") +
#   annotate("text", x = as.Date("2020-06-01"), y = max(decomp_df$observed)*0.8,
#            label = "COVID-19 Impact phase\n(2020–2023) ", size = 2.8, color = "#333333")

## 趋势分量
p_trend <- ggplot(decomp_df, aes(x = date)) +
  geom_line(aes(y = trend), color = "#D62728", linewidth = 0.8) +
  labs(x = NULL, y = "Trend Component", title = "Long-term Trend") +
  idp_theme+
  scale_x_date(
    breaks = seq(ymd("2005-01-01"), ymd("2023-01-01"), by = "2 years"),
    date_labels = "%Y",
    expand = expansion(add = c(30, 60))
  ) 

## 季节分量
p_seasonal <- ggplot(decomp_df, aes(x = date)) +
  geom_col(aes(y = seasonal), fill = "#2B8CBE", width = 20, alpha = 0.7) +
  labs(x = NULL, y = "Seasonal Variation", title = "Seasonal Pattern") +
  idp_theme +
  scale_x_date(
    breaks = seq(ymd("2005-01-01"), ymd("2023-01-01"), by = "2 years"),
    date_labels = "%Y",
    expand = expansion(add = c(30, 60))
  ) 
  # geom_rect(xmin = as.Date("2020-03-01"), xmax = as.Date("2020-08-01"),
  #           ymin = -Inf, ymax = Inf, fill = "#FFEDA0", alpha = 0.15) +
  # annotate("text", x = as.Date("2020-05-15"), y = max(decomp_df$seasonal)*0.8,
  #          label = "High-risk season\n(Mar-Aug)", size = 2.8, color = "#333333")

## 残差分量
p_residual <- ggplot(decomp_df, aes(x = date)) +
  geom_line(aes(y = remainder), color = "#2B8CBE", linewidth = 0.4) +
  geom_hline(yintercept = 0, color = "#D62728", linetype = "dashed") +
  labs(x = NULL, y = "Residuals", title = "Random Residuals") +
  idp_theme +
  geom_point(data = filter(decomp_df, abs(remainder) > 10),
             aes(y = remainder), color = "#E31A1C", size = 1.5)+
  scale_x_date(
    breaks = seq(ymd("2005-01-01"), ymd("2023-01-01"), by = "2 years"),
    date_labels = "%Y",
    expand = expansion(add = c(30, 60))
  ) 

# 组合图形 ----------------------------------------------------------------
final_plot <- (p_original / p_trend / p_seasonal / p_residual) + 
  plot_layout(heights = c(1.2, 1, 1, 1)) +
  plot_annotation(tag_levels = 'A', tag_suffix = ')') &
  theme(plot.tag = element_text(face = "bold", size = 10),
        plot.margin = unit(c(5, 10, 5, 5), "mm"))
final_plot

# 保存输出 ----------------------------------------------------------------
ggsave("./主图/时间序列分解图.tiff", 
       plot = final_plot,
       width = 21, 
       height = 28,
       units = "cm",
       dpi = 600,
       compression = "lzw")

# 导出分解数据 ------------------------------------------------------------
write.xlsx(decomp_df, "./output/时间序列分解成分.xlsx")
